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NEAR-OPTIMAL PERFECTLY MATCHED LAYERS 
FOR INDEFINITE HELMHOLTZ PROBLEMS 

VLADIMIR DRUSKIN*, STEFAN GUTTELt, AND LEONID KNIZHNERMAN* 

Abstract. A new construction of an absorbing boundary condition for indefinite Helmholtz 
problems on unbounded domains is presented. This construction is based on a near-best uniform 
rational interpolant of the inverse square root function on the union of a negative and positive real 
interval, designed with the help of a classical result by Zolotarev. Using Krein’s interpretation of 
a Stieltjes continued fraction, this interpolant can be converted into a three-term finite difference 
discretization of a perfectly matched layer (PML) which converges exponentially fast in the number 
of grid points. The convergence rate is asymptotically optimal for both propagative and evanescent 
wave modes. Several numerical experiments and illustrations are included. 

Key words. Helmholtz equation, Neumann-to-Dirichlet map, perfectly matched layer, rational 
approximation, Zolotarev problem, continued fraction 
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1. Introduction. An important task in science and engineering is the numerical 
solution of a partial differential equation (PDE) on an unbounded spatial domain. 
Unbounded spatial domains need to be truncated for computational purposes and 
this turns out to be particularly difficult when the PDE models wave-like phenomena. 
In this case the solution may not decay rapidly towards the truncation boundary and 
artihcial reflections and resonances may pollute the numerical solution. The prototype 
of such notorious PDEs is the Helmholtz equation, which models the propagation of 
electromagnetic or acoustic fields from a source with a single frequency k > 0 

c^Au + k^u = 0. (1.1) 


As a motivating example from geophysics, this equation may be posed on an un¬ 
bounded half-space corresponding to the Earth’s subsurface and the variable wave 
speed c may be caused by variations in the sedimental composition, see Figure 0 
In seismic exploration a pressure wave signal of frequency k is emitted by an acoustic 
transmitter placed on the Earth’s surface or below, travels through the underground, 
and is then logged by receivers. From these measurements geophysicists try to infer 
variations in the wave speed c which then allows them to draw conclusions about 
the subsurface composition. Clearly, the spatial domain for this problem needs to be 
truncated and there are various ways for achieving this, with a very popular approach 
being known as perfectly matched layer (PML, see [7 10 13 ). 

A perfectly matched layer can be seen as a localized modification of the spatial 
discretization scheme to absorb the waves exiting the computational domain. In a 
finite difference framework such layers typically lead to variable complex-valued step 
sizes, which is why this approach is sometimes also referred to as complex coordinate 
stretching. The aim of an efficient PML is to achieve a strong absorption effect by 
adding only a few number of layers. The aim of this work is to extend a modern 
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Figure 1.1. A typical setup in seismic geophysical exploration where a source emits pressure 
waves into the Earth’s subsurface which are then logged at (multiple) receivers. The wave propagation 
is modelled by the Helmholtz equation The varying shades of gray in the Earth’s underground 

indicate different values of the wave speed c = c{x,y). The perfectly matched layers developed in 
this paper allow for a variation in c tangential to the boundary of the computational domain. We 
will return to this example in section\6.S\ 


finite-difference construction of perfectly matched layers which are near-optimal for 
indefinite Helmholtz problems, that is, they achieve near-best possible absorption for 
a given number of layers. The number of required layers is critical in particular for 
large-scale simulations of three-dimensional exterior problems. A variety of such prob¬ 
lems arise, for example, in oil and gas exploration, and near-optimal grids are part of 
almost all electromagnetic simulators used at Schlumberger [I 15,57. Other appli¬ 


cations of effective discretizations of exterior domains include homogenization theory, 
photonic crystals, energy-driven pattern formation, and the modelling of biologic cell 
communication (see, e.g., 41-4^). 


1.1. Outline of this work. We will now give a short overview of this work 
and explain the structure of the paper. Let us start by considering a prototype of a 
differential equation on an unbounded domain, the two-point boundary value problem 


d'^ , d 

u = Au, 


dx^ 


dx '“-0 


= -b, 


u 


a:=+oo 


= 0 , 


( 1 . 2 ) 


where A G 


■^NxN 


is nonsingular and {b,u{x)} C C^. If A is a discretization of a 


differential operator on some spatial domain H C then (1.21 is a semidiscretization 
of an {£ + l)-dimensional partial differential equation on [0, -foo) x H. Assuming that 
problem (1.2) is well posed (which may require some additional conditions like, e.g., 
the limiting absorption principle discussed below), its exact solution can be given in 
terms of matrix functions as u{x) = exp(—In particular, at a; = 0 the 
solution is given as 


u{0)=F{A)b, F{z) = z 


= .- 1/2 


(1.3) 


The function F{z) is often referred to as the impedance function (also known as Weyl 
function), and it completely characterizes the reaction of the unbounded domain to an 
external force [^. The relation (1.3) allows for the exact conversion of the Neumann 


data — b at the boundary x = 0 into the Dirichlet data it(0), without the need for 


solving (1.2) on its unbounded domain. This is why T(A) is often referred to as the 


Neumann-to-Dirichlet (NtD) operator. 
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When solving wave scattering problems one typically deals with a discretization 
matrix A of the negative shifted Laplacian — c^A — on fl C Under the assump¬ 


tion that c does not depend on x, problem (1.2) is a semidiscretization of the indefinite 
Helmholtz equation on the domain [0, -boo) x fl. In this case the matrix A is 
(similar to a matrix) of the form 


A = L-k^I, 


(1.4) 

is the identity matrix, 
and > 0 is not in the spectrum of L. For a solution of (1.2) to be unique we impose 


where L C 


^NxN 


is Hermitian positive definite, I S 


tiNxN 


the limiting absorption principle (see, e.g., 50 ). This means that for a real number 
k we define it as a limit of solutions it(^+®‘^) of (1.21 with wave numbers k + ie (e > 0) 
instead of fc, i.e.. 


It = lim 

e\,0 


(1.5) 


This uniquely defines the value F{A) = notwithstanding that some eigenvalues 

of A may lie on the standard branch cut of F{z). 

We will now outline our construction in the following sections, which combines 
ideas of the eminent mathematicians Y. I. Zolotarev (1847-1878), T. J. Stieltjes (1856- 
1894), and M. G. Krein (1907-1989). The main aim in section i is to approximate 
F{z) by a rational interpolant Rn{z) of type (n — 1, n), so that Rn{A) can be seen as 
an approximate NtD operator, mapping the Neumann data —6 to the Dirichlet data 
Rn{A)b. Clearly, the weighted 2-norm approximation error of this map is 


|/i,.(4)6 - F(A)(,|| = 


\ 


N 


^|i?„(A,-P)-F(A,-F)P|6, 


i=i 


where bj = v*b and are the eigenpairs of L with ||vj|| = 1. We have Ai < 

< Aat and thus arrive at the problem of scalar rational approximation of F(z) 
on the union of a positive and a negative real interval. Our rational interpolant 
Rn{z) is obtained by combining two optimal Zolotarev interpolants constructed for 
the two intervals separately. For illustration purposes we have graphed the relative 
error of such a function in Figure 1.2 In addition to the explicit construction of such 


approximants, section also contains a novel detailed convergence analysis, with the 
more technical proofs given in the appendix. 

In section]^ we will show that the rational function i?„(z) can be converted into 
an equivalent three-term finite difference scheme on a nonuniform grid with n points. 
This is achieved by formally rewriting Rn{z) as a Stieltjes continued fraction and 
using Krein’s interpretation of that fraction as a finite-difference scheme. However, 
due to the non-Stieltjes nature of i?„(z) (its poles may lie on a curve in complex 
plane, as shown in Figure [2T] ) the continued fraction coefhcients can also be complex, 
which results in a finite difference scheme with complex-valued grid steps. This scheme 
allows for the simple and efficient computation of an NtD map and the construction of 
an absorbing boundary layer for indefinite Helmholtz problems. The near-optimality 
of Rn{z) implies that the number of required grid points is close to smallest possible. 
A summary of an algorithm for computing this grid is given in section 

Section is devoted to the adaptation of our PML construction to a second-order 
finite difference framework. In particular, in section [5T| we extend our optimal ratio¬ 
nal approximation approach to the infinite lattice problem. Our analysis carries over 
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Figure 1.2. Relative error \RrL(z)/F(z) — 1| of a rational approximant Rn(z) for F{z) = z~^/'^ 
on [—le3, —1] U [1, le4]. We have adopted a special plotting type for simultaneously visualizing large 
intervals on the negative and positive real semiaxes in logarithmic scales, with the gray linear region 
in the middle gluing the two intervals together. The rational function Rn{z) is of type (n — l,n), 
n = 9, and it has been constructed by combining two Zolotarev interpolants with mi = 8 and m 2 = 10 
interpolation nodes for the negative and positive intervals, respectively. Visually the solution of our 
complex rational approximation problem behaves similarly to the max-norm optimal errors of the real 
problems, i.e., it shows “equal ripples” on the targeted i nter vals (although the Chebyshev alternation 
theory Ch. II] is not applicable in the complex case [52]). 


to this problem, thereby providing a novel theoretical justification for the exponential 
error reduction in our perfectly matched layer as the number of grid points increases. 

Finally, in section we demonstrate the high accuracy and exponential conver¬ 
gence of our perfectly matched layer with several numerical examples. 


1.2. Review of related work. It was already shown in 18 34 that a rational 
approximant Rniz) of type (n — l,n) for the function F(z) can be converted into 
an equivalent three-term hnite difference scheme on a special nonuniform grid with 
n points, mapping the Neumann data —b to the Dirichlet data Rn{A)b. In these 
papers the authors were mainly concerned with a special instance of (1.2) where A 
corresponds to a discretization of the negative Laplacian —A, in which case A is 
a real symmetric positive definite matrix. The error of the approximate Neumann- 
to-Dirichlet (NtD) map is then bounded by the maximum of \Rn{z) — F(^z)\ on the 
positive spectral interval of A. Approximation theory allows for the construction 
of exponentially convergent rational functions Rn{z) with a convergence rate weakly 
dependent on the condition number of A, thus producing a three-term finite difference 
scheme with a so-called optimal grid (also known as finite-difference Gaussian rule or 
spectrally matched grid). The connection of Rniz) and this grid is inspired by Krein’s 

It was shown in 


mechanical interpretation of a Stieltjes continued fraction 36 


that the same grids produce exponentially convergent NtD maps even for problems 
arising from the semidiscretization of anisotropic elliptic PDFs and systems with 
mixed second-order terms, i.e., when the second-order ODE system in (1.2) is modified 
by adding a first-order term. 

It should be noted that the positive and negative eigenmodes of A correspond 
to so-called evanescent and propagative solutions riiA)uix) and r]i—A)u{x), respec¬ 
tively, with 77(5) denoting the Heaviside step function. The evanescent modes, i.e., the 
nonzero eigenmodes in the spectral decomposition of 77 ( 2 !) M(a;), decay exponentially 
as X increases (hence the name). Therefore a simple, though possibly not the most 
efficient, way to absorb them is to truncate the domain at some (sometimes quite 
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significant) distance from the targeted area of interest, and then to deal with the 
propagative modes alone. On the other hand, the norm ||r 7 (—A)u(a:)|| does actually 
not depend on x, so simple boundary truncation will not be effective for absorbing 
propagative modes. 

Enquist and Majda computed Rn{z) as a Fade ap- 


In their seminal paper 21 


proximant of F(z) at some real negative point and then evaluated it via continued 
fraction-type recursions. This approach yielded exponential convergence on the neg¬ 
ative real semiaxis, however, with the rate quickly deteriorating towards the origin. 
Another celebrated approach for absorbing propagative modes is called complex scal¬ 
ing and was originally introduced in for molecular physics calculations. It is also 


known as perfectly matched layer (PML), a term coined in the influential work 10 


where it was independently rediscovered and adapted for time-domain wave propa¬ 
gation. We will use the latter term because it seems to be more established in the 
wave propagation literature. The well-posedness of the PML formulation was studied 
in [^[^. The essence of the PML approach is a complex coordinate transformation 
which changes purely imaginary exponentials of propagative modes to complex decay¬ 
ing ones, thus, in principle, allowing reflectionless domain truncation mu- However, 
coarse PML discretizations introduce undesirable numerical reflections which decay 
rather slowly with the grid size in case of low-order discretization schemes. This 
problem was partially circumvented in for the solution of time-domain wave prob¬ 
lems, where the optimal gridding approach was extended to PML discretizations. By 
choosing an appropriate purely imaginary grid this approach allowed for the construc¬ 
tion of all possible rational interpolants Rn{z) for F{z) on a real negative interval, 
including the Fade approximants constructed in [21] , and preferably the best uniform 
approximants targeting the spectral support of the expected solution. See also 39 
and 19 for adaptations of the optimal gridding approach to the hyperbolic elastic¬ 
ity system and the Helmholtz equation, respectively. A non-optimal PML layer for 
absorbing both evanescent and propagative modes in dispersive wave equations has 
been proposed in 56 . However, the problem of designing discrete PMLs which are 


optimal for both wave modes remained open. 

Our construction in section]^ is inspired by a “trick” originally used by Zolotarev 
and Newmann, writing the relative approximation error Rn(z)/F{z) — 1 in terms of 

i(s) is a polynomial of degree m = 2n, = z. This trick 


Hm{s)/Hm{-s), where H, 
was rediscovered in 


28,^, where F[ra{s)/F[m{—s) was identified with the numeri¬ 


cal reflection coefficient, and a continued-fraction absorbing condition was explicitly 
constructed in terms of the roots of Hm{s) and introduced in the PDE discretization 
via a so-called trapezoid finite element method. However, these important papers fell 
short of introducing optimal approximants. In addition to the construction of these 
approximants, section[^also contains a novel detailed convergence analysis. To make 
our paper more pleasant to read we have decided to present the technical proofs in 
an appendix. 

In an unfinished report 17 , the authors suggested to split i7m(s) into the product 
of polynomials with real and imaginary roots, thus decoupling the approximation 
problems on the positive and negative intervals. It was then suggested to apply 
conventional optimal rational approximants on each of the two intervals, and the 
resulting error was only determined by the largest error of these two approximants. 
A drawback of such an approach is that it requires the splitting of the PML grid into 
two subdomains with nonlocal finite difference stencils at the conjugation interfaces. 
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2. Construction of a near-optimal approximant on two intervals. The 

function is commonly defined in the complex plane C with the slit (— 00 , 0] C K. 

However, in our application we need an analytic continuation F[z) of 2 “^/^ from the 
positive real semiaxis M+ = {cc S K | a: > 0 } to —K+ in accordance with the limiting 
absorption principle (1.5), i.e., attaining the values 

F{z) = -i{-z)-^/^ 


for 2 G -I 




and the principal value of the square root for 2 G M+. We will therefore assume in 
the following that F{z) is defined in C with the branch cut in the lower half-plane. 

Following [^, we now construct a rational interpolant Rn{z) of type in — l,n) 
to F{z) on the union K of two real intervals 

K = KiUK2, iFi = [ai,6i], K 2 = [a2,b2], 
oi < < 0 < 02 < 62 , 

using solutions of a classical Zolotarev problem on each of the two intervals. In view 
of the definition (1.4), these intervals will correspond to the spectral subintervals 
[Ai — Xig — /c^] and — k'^, Xn — k"^] (or their estimates), respectively, where 

Ai < • • • < Ai„ < < Aig+i < • • • < Xn- 

Separating the odd and even parts of a polynomial Hm of degree m = 2n, we 
define polynomials Pn-i and Qn of degrees < u — 1 and n, respectively, such that 


Hmis) = -sP„-i(s^) + Q„(s^). 


The rational function 


Rniz) = 


Pn-ljz) 

Quiz) 


( 2 . 1 ) 


( 2 . 2 ) 


will be considered as an approximant for F(z) on K. We have 

rt I 2\ _ ^ Pn-l{s ) _ H^[ — s) — 

Q„(s 2 ) - H^i-s) + H^{s) ^ ^ 

and thereby obtain an expression of the relative averaged approximation error as 


\F{s^) - R„is^)\ ^ 

Hm{s) 

|F(s 2 ) + i?„(s 2 )| “ 

Hm{-S) 


Following [17[ Section 2], we can split the approximation problem on K into two 
independent problems on Ki and K 2 . 

Lemma 2.1. Let mi and m 2 be positive integers such that m = mi -I-m 2 , and let 
Hmi CLnd Hm 2 be polynomials of degrees mi and m 2 with roots on F{Ki) and F[K 2 ), 
respectively. Define 




Then 


H^{s) 

= max 

^mi (■s) 

H^-s) 

( 5 ) 

sGF(Ai) 


max 

''eF{K^) 
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and 


max 

seF{K2) 


Hmjs) 

Hmi-S) 


max 

s£F{K2) 


Hm2 (s) 
Hm2 {~s) 


Proof. This lemma immediately follows from the equalities 


dp mi i^) 
dim I ( 


1 if s e f{K2), 


and reciprocally 


dim2 (s) 
ddm2 ( s) 


1 if s e F{Ki). 


Let us consider a single real interval [c, d] with 0 < c < d, and the problem of 
finding a real monic polynomial Zm‘^'^ of degree m > 1 (denoted as Zm'^'^ S 'Pm,reed) 
which attains the minimum in the Zolotarev problem 


= min max 

ZGPm.roal C<S<d 


Zis) 


Zi-s) 


(2.4) 


It is known from 


4^ 58 that this minimizer Zm'^'^ exists uniquely, that its roots 
(j = 1,..., m) are locaied in (c, d), and that they are expressible in terms of elliptic 
integrals. More details are given in the appendix, in particular, formula (A.l). 

We choose positive integers m-i and m 2 and introduce the polynomial 

Hm{s) = (2.5) 

of degree m = mi + m 2 . From Lemma |2.1| we obtain the following result. 


Proposition 2.2. The polynomial FIm{s) defined in (2.5) satisfies 
Hmis) 


max 

s£F{K) 


Hmi-S) 




It is well known that the classical Zolotarev functions in (2.4) converge expo¬ 


nentially. Let us denote by the Cauchy-Hadamard convergence rate of Zm^'^ 


i.e.. 


= lim S = c/d. 


An exact expression of p^^^ in terms of elliptic integrals is given in (A.3). For small 
interval ratios S one can derive a simple approximate expression 




exp 




in terms of elementary functions 34, Appendix A]. This expression shows the weak 
dependence of the Cauchy-Hadamard convergence rate on the interval ratio S. 
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In view of Proposition |2.2| mi and m 2 should be chosen to balance the errors of 
both Zolotarev functions. One way of achieving this is by setting 


Pi =p(v^), 


P 2 = p(v'W6l) 


and 


mi = m 


log P 2 


log Pi + log P 2 


+ 0, m2=m —mi, |0| < 1/2, 


( 2 . 6 ) 


(2.7) 


where 0 is chosen to round (mlogp 2 )/(log/ 9 i +logp 2 ) to the nearest integer. We are 
now in the position to formulate a near-optimality result for the obtained approximant. 

Theorem 2.3. Let us denote 


p = exp 


log Pi logP 2 
log Pi + log P 2 


( 2 . 8 ) 


Let the polynomial Hm be defined by (2.5), the polynomials Pn-i and Qn defined 


by (2.1), the rational fraction i?„ defined by (2.2), and m = 2n. Further let the 

(2.9) 


conditions (2.7) and 

2max|p7^/^,p^^/^|p'" < 1 
he satisfied. Then the upper relative error bound 


max 

zgk 


R,fiz) 

F{z) 


- 1 


< 


4max|pi ^/^p2 
1-2 max p" 


( 2 . 10 ) 


holds. On the other hand, if P and Q ^ 0 are arbitrary polynomials of degrees < n — 1 
and < n, respectively, then R = P/Q satisfies the lower error bound 


max 

zgk 


R{z) 


F{z) 


- 1 


> 


2p" 


1 -hp” 


( 2 . 11 ) 


This theorem, whose proof is given in the appendix, implies that the upper er¬ 
ror bound for our Zolotarev approximant Rn{z) and the lower bound for the best 
possible approximant have the same Cauchy-Hadamard convergence rate p, i.e., our 
approximant is asymptotically optimal in the Cauchy-Hadamard sense. As is also 


demonstrated by the following numerical example (and the corresponding Table 2.1) 


the Zolotarev approximant can be worse than the best possible approximant only by 
a moderate factor. We should point out that, unlike their real counterparts, complex 
max-norm optimal rational approximation problems are generally not convex and may 
have non-unique solutions 52 . It therefore seems unlikely that the near-optimality 


result of Theorem |2.3| can be improved significantly. 

Example 2.1. Let us, as in Figure consider the problem of approximating 
F(z) = z~^l'^ by a rational function Rn{z) of type (n — l,n) on the union of two 
intervals K = [ai, 6 i] U [ 02 , 62 ] = [—le3,—1] U [1, led]. Using the exact formula (A.3| 
we calculate 


Pi « 0.361, p2 ~ 0.439, p « 0.634. 
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Tab le 2 .1 

Lower and upper error bounds of Theorem \2.3\ and actual errors max^gx \Rn{z)/F{z) — 1|, 
F{z) = z~^/'^, for various values of m = 2n. The set K is chosen as K = [—le3, —1] U [1, le4]. 


m 

mi 

m 2 

bound (2.11 

1 

relative error 

bound ( 

2.10 

6 

3 

3 

1.22e- 01 

3.42e - 01 

5.52e- 

-01 

12 

5 

7 

8.41e- 03 

2.47e - 02 

2.85e- 

-02 

18 

8 

10 

5.49e - 04 

1.15e- 03 

1.83e- 

-03 

24 

11 

13 

3.57e- 05 

8.95e - 05 

1.19e- 

-04 

30 

13 

17 

2.32e - 06 

7.01e - 06 

7.72e - 

-06 

36 

16 

20 

1.51e- 07 

3.29e - 07 

5.02e- 

-07 

42 

19 

23 

9.79e - 09 

2.37e - 08 

3.26e - 

-08 

48 

21 

27 

6.36e- 10 

2.01e- 09 

2.12e- 

-09 

54 

24 

30 

4.13e- 11 

9.43e - 11 

1.38e- 

- 10 

60 

27 

33 

2.69e- 12 

6.28e- 12 

8.94e - 

- 12 


In Table \2A\ we list the error bounds of Theorem 2.3 for various values of m = 2n 
together with the actual approximation error. The calculations confirm the bounds 
and show that they are roughly of the same order, i.e., our approximants Rn{z) have 
relative errors of the same order as the best possible approximants. 


The logarithmic surface plot in Figure 2.1 shows the relative error \Rn[z)/F{z) — V\ 
for the case n = 9 (the same as in Figure 1.2). Note how the poles align on a 


curve in the lower-left quadrant of the complex plane. We speculate that this curve 
asymptotically (as n —>■ oo) approximates the shifted branch cut C of the analytic 
continuation of F{z) into the lower half-plane, and that C possesses the so-called S- 
property (“symmetry property”, see ^6 with respect to K. This would imply 

that the equilibrium charge of the condenser {K, C) has a logarithmic potential which 
is (constant and) minimal on K over all “attainable” branch cuts. Our experiments 
also suggested that the curve C coincides exactly with the negative imaginary semiaxis 
in the case of symmetric intervals Ki = —K 2 , and that it approaches the real positive 
or negative semiaxis for large or small ratios mi/m 2 , respectively. 


A remarkable feature in Figure 2.1 is that the relative error \Rn{z)/F(z) — 1 | 
stays uniformly small “above” the set K, i.e., for complex numbers z with positive 
imaginary part and real part in K. We will return to this observation in section\67^ 


3. Finite difference grids from rational approximants. We now explain 
how a rational function i?„(z) « F{z) can be transformed into an equivalent staggered 


finite difference grid for (1.2). Assume that we are given primal grid points and steps 


0 — xq, Xi, ..., x^i, hj — Xj Xj—i, 


and dual grid points and steps 


0 — Xq, Xi, ..., Xji, hj—i — Xj Xj—i, 


with j = 1,..., n in both cases. Denote by Uq, Ui,..., Un approximations to the 


solution u{x) of (1.2) at the primal grid points xo,Xi,... ,Xn. Let the first-order 


finite differences {uj — Uj-i)/hj be located at the dual points % (j = 1,... ,n). We 
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real(z) 


Figure 2.1. Relative error \Rn(z)/F(z) — 1| of a Zolotarev approximant Rniz) for K = 
[—le3, — 1] U [1, le4] and n = 9 shown as a logarithmic surface plot over a region in the complex 
plane. The imaginary axis is plotted in reversed direction for a better panoramic view. 


assume that the following finite difference relations 


/ Ui- Up 
hp \ hi 
1 f Mj + l — Uj Uj — Uj-i 


y+i 


+ 6 ) = Aup, 

= Auj, j = 


(3.1) 

(3.2) 


are satisfied with the convention that = 0. It can be verified by back-substitution 
that the value tip specified by these recursive relations can be written as 


Up = Rn{A)b, 


where Rn{z) is a rational function of type (n — l,n). By construction, —Rn{A)~^ is 
the Schur complement of the submatrix with positive indices of the system (3.1 )-(3.2). 
Written as a finite-length Stieltjes continued fraction (S-fractior0 this function takes 
the form 


Rn{z) = 


hpz + 


hi 


hiz - 


hn-l + 


hn-lZ+ — 
hn. 


(3.3) 


Recalling from above that the exact solution of (1.2) satisfies tt(0) = A we are 

apparently left with the problem of determining Rn{z) such that Rn{A)b « 


^We now allow for complex-valued hj—i^hj {j = 1,. .. ,n) in (3.3i, which is different from the 
classical definition of S-fractions with real positive parameters. 
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optimally in some sense. The conversion of Neumann data —ft to Dirichlet data tt(0) 
can now be realized by solving a finite difference relation on a grid generated from 


quantities hj-i and hj (j = 1,... ,n) in (3.3). 


33 


The connection between the S-fraction (3.3) and the finite difference problem 
(3.1)-(3.2| is due to Mark Krein (see, e.g., |36|). He viewed the problem (3.1|-(3.2| 
as a so-called Stieltjes string, which is a string of point masses /ij-i and weightless 
stiffnesses hj {j = 1,... ,n), both real positive. There is a one-to-one correspondence 
between the set of Stieltjes strings and Stieltjes spectral functions Rn{z), which are 
rational functions of type (n—1, n) having n non-coinciding real negative poles and real 
positive residues. For this case, the S-fraction parameters hj_i and hj {j = 1,... ,n) 
can be computed via 2n steps of the Euclidean polynomial division algorithm (see, 
which can be stably executed with the help of the reorthogonalized Lanczos 
The optimal rational approximation of F(z) on a positive real interval 
is a Stieltjes problem 34 , hence the generated grid steps are real positive. The 
approximation problem on a single negative interval can be solved by using i?„(—z), 
where Rn{z) is the approximation on the symmetrically reflected positive interval. 
This reflection rotates the grid steps hj-i and hj (j = 1,..., n) by an angle of 7r/2 in 
C, i.e., it makes the grid steps purely imaginary. Generally, the problem of optimal 
approximation on the union of a positive and a negative interval leads to non-Stieltjes 
rational functions Rn{z) of type (n — l,n). Assuming absence of breakdowns (which 
are unlikely but can not be definitely excluded), the transformation to the non-Stieltjes 
rational function (3.3) can still be carried out via the complex 2n-step Euclidean 
algorithm. 


e.g. 

algorithm 18 


We used the bi-Lanczos extension of the Lanczos-based algorithm [18| 
which, according to our experience, always produced meaningful results. 

Example 3.1. We begin with reproducing a real optimal grid from (3.3) generated 
for a real positive interval K = [l,le4], see Figure 3.1 (left). Similar results were 
reported in 34 . We can consider this example as a degenerate case of the two-interval 
problem with mi = 0 and m 2 = 10. The plot shows “alternation” of the primal and 
dual grid points and monotonically growing steps. The grid looks like an equidistant 


grid stretched by a rather smooth transform. R was shown in 34 that for large n and 


small interval ratios such transforms are asymptotically close to the exponential. 


In Figure 3.1 (right) we plot the complex finite difference grid points obtained 
from the continued fraction (3.3) in the case when K = [—le3, —1] U [l,le4] and 
mi = 8 and m 2 = 10. We notice the “alternation” of the primary and dual points 
on some “curve”, which is an intuitive evidence of a good quality of the grid, i.e., we 
can speculate that the finite difference solution approximates the exact solution with 
second-order accuracy on that curve. This curve can be interpreted as the complex 
PML transform of the real positive axis in accordance with [7 13 . 


In summary, we observe that the finite-difference operators on grids obtained 
from (3.3) approximate the second-order derivative operator on curves in the com¬ 
plex plane. This can be viewed as a complex extension of Krein’s results on the 
convergence of the Stieltjes discrete string with impedance Rn{z) to its continuous 
counterpart with impedance F{z) when —>■ E on IR+ 36 . Besides internal beauty. 


this phenomenon may have useful consequences. Eor example, it lets us hope that 
pseudospectral estimates and stability results for continuous PMLs and damped ID 


differential operators kyD 16 remain valid for (3.1)-(3.2) with the optimal grid. 


4. Summary of the algorithm. In the following we provide a step-by-step 
description for computing the grid steps hj-i and hj {j = 1,... ,n) in (3.3). 
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□ primal grid points 
O dual grid points 
















10 ^ 10 ^ 10 ‘ 10 “ 10 ' 



Figure 3.1. Grid points generated from quantities in the continued fraction ( |3.3[ | . Left: In this 
single-interval case the set K is chosen as K = [l,le4] with m 2 = 10 (and mi = 0^. Right: The 
set K is chosen as K = [— le3, —1] U [l,le4] with mi = 8 and m 2 = 10. The gray ^‘continuous” 
curve has been obtained by connecting the grid points generated with the parameters mi = 27 and 
m 2 = 33, and we conjecture that the grid points align on a limit curve as m ^ 00 . 


1. It follows from (2.3) that the numbers 




W-bi,\/-<^i) 
j 


) , j = and ^ j = l,...,TO 2 , 

are the interpolation nodes for Rn{z) as an interpolant of F{z). Knowing in¬ 
terpolation nodes and function values, we compute the coefficients of P„_i(z) 
and Qn {z) by means of solving the corresponding system of linear algebraic 
equations in high-precision arithmetic. 

2. The poles of the interpolant, i.e., the roots of Quiz), can be computed as the 
eigenvalues of an associated companion matrix, see [M 
solve this eigenvalue problem we use the quasi-versior 


Subsection 7.4.6], To 
^ of the QR transfor¬ 
mation method 45, § 11.6] and then, if necessary, correct the roots by means 
of a combination of Laguerre’s 45 ,§ 9.5] and Newton’s 37 method. 


Knowing the poles of Rn{z), the corresponding residues are computed. 
Finally, the grid steps hj-i and hj (j = 1,..., n) are computed using the re¬ 
cursion formulas [18[ (3.4)], with the underlying analogue of an inverse eigen¬ 
value problem for a symmetric tridiagonal matrix (see [18[ subsection 3.1, 
item 3°], 46 theorem 7.2.1]) being solved by a quasi-Lanczos process 


14 


Ch. 6] with quasi-reorthogonalization. Here we used the well-known connec¬ 
tion between the Lanczos and Euclidean algorithms (see, e.g., 33 ). 


5. Adaptation to a second-order finite difference framework. 

5.1. Approximation of the discrete impedance function. So far we have 
considered the function E(z) = which arises when solving the boundary-value 


problem (1.2) for x S [0, -foo). When this problem is seen as an infinite extension of 
some interior computational domain, the exponential convergence of the interpolant 
Rn{z) is consistent with a high-order (or even spectral) discretization of the operator 
acting in this computational domain. 

However, it is also possible to compute the NtD map of a discretized version 


of (1.2) on a uniform infinite grid via rational approximation of a slightly modified 


we formally use in the complex case the formulas intended for the real case. 
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function Fh{z) to be determined below. This function will lead to a three-term finite 
difference scheme which is appropriate for being combined with a standard second- 
order finite difference discretization in the interior computational domain, because it 
allows for the elimination of spurious reflections from the PML boundary due to the 
error of the interior discretization. 

Given a fixed step size h > 0, let us consider the problem (3.1)-(3.2) on the 


infinite equidistant grid with /ig = 0.5h and hj = hj = h for j = 1,..., oo. We will 
determine a function Fh{z) such that 

tiQ = Fh{A)b 

via a well-known approach widely used in the representation of irrational numbers 


51 


via continued fractions (see, e.g., 20 section 9]). This approach was already applied 
to the infinite lattice problem: the infinite-length S-fraction representation of 


Fh analogous to (3.3) is 


Rh{z) = 


0.5hz 


h + 


hz ■ 


(for a proof of convergence we refer to 49 
continued fraction 


or 


35, theorem 4.58]). The remainder 


S{z) = 


hz + 


hz - 


evidently satisfies the equation 


S{z) = - - - 

h —r 

hz + S{z) 


or equivalently S(z)‘^ + hzS{z) — z = 0. Since 0.5hz + S{z) = i?„(z) ^ = Fh{z) we 
have arrived at the quadratic equation 


Fhizf 


1 

z+{0.5hz)^' 


We choose the root which converges to the exact impedance F{z) as /i —>■ 0, i.e.. 


Fh{z) 


1 

\/ z + (O.bhz)^ 


(5.1) 


This function, which we will refer to as the discrete impedance function, approximates 
with second-order accuracy the exact impedance at the boundary, so being centered, 
the resulting finite difference scheme is of second order globally. 
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Analogously to what we had achieved with (1.3) for continuous x, the relation 


(5.1) allows us to convert the Neumann data — b at a; = 0 into the Dirichlet data ttg 


without actually solving the infinite lattice problem. 

7 2 

For a given b > 0 let us define cr = ^. The invertible linear fractional change of 
variables 


CTZ + 1 


(5.2) 


translates the union of a negative and a positive segment K = [oi, bi] U [ 02 , ^ 2 ] again 
into the union of a negative and a positive segment. Let us assum^that < ai. 

Let Pn-ilQn denote the rational approximant of theorem |2.3| for the image of K 


under transformation (5.2). Then 


Pn-i{w) 

w ■ ——-—^-1 




(TZ+l Qni-^)icrz+l) 


Qn{w) 


- 1 


is small on K, the numerator and the denominator 


Fhiz)■ 


az + 1 
^n-i(yy+T 




Q"(o-z+l) 

)(ctz + 1)"-i 


-1 


Qri(-^^fj)(0'Z + 1 )’’ 


- 1 




n—1 


az + 1 


(az + 1) 


n—1 




az + 1 


(az + 1)’^ 


being polynomials of degrees < n — 1 and < n, respectively. We have thereby es¬ 
tablished a direct relation between the errors of the rational interpolants for P(z) 
and Phiz) on transformed compact sets, respectively, with the interpolation nodes 
being transformed accordingly. This allows us to conclude that we obtain identical 
convergence rates for both interpolation processes. In particular. Theorem |2.3| holds 
with F{z) being replaced by Fh(z). 

We would like to mention that a rational approximation-based absorbing bound¬ 


ary condition for the infinite lattice was suggested in 51 and combined with a trape¬ 


zoidal finite element approach in 29 . However, that approach required a modification 


of the Helmholtz equation by a higher-order term. On the contrary, in our framework 
the discreteness can be incorporated simply by adjusting the PML grids. Visually 
these grids look very similar to the ones shown in Figure |3.1[ i.e., we can specu¬ 
late again that they approximate the exact solution u(x) of ( |1.2[ ) with second-order 
accuracy on some modified x-curve in the complex plane. 

5.2. Matching interior and exterior discretizations via a single grid. 

Let us consider the second-order infinite equidistant finite difference problem 


%+i 


h 


'^ 3-1 

h 


AUj — Qj, 


J = 


with boundary conditions 


= 0 , 


lim Uj = 0, 

J-fOO 


,-1,0,1,...,00 (5.3) 


(5.4) 


^As discussed earlier, the parameter ai should be set to a lower bound of A’s spectral interval, 
in which case the condition —cr~^ < ai corresponds to the Nyquist sampling criterion of two grid 
points per wave length. This assumption should be met by any reasonable discretization scheme. 
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assuming qj = 0 for j > 0. Problem (5.3) can be split equivalently into an interior 
finite-dimensional system 


1 ( Mj + l — Itj Ml — 


- 


J = 




(5.5) 


1 


0.5/i \ h 

and an exterior infinite system 

1 f Ui — Uq 


^ ^ ) - Auo = 0, 


0.5h \ h 

1 f — Uj Uj — Uj_i 
h h 


b ) — Auq = 0, 


(5.6) 


— Aui = 0 , 


j = l,...,oo, 


both systems being coupled via a vector variable 

Problem (5.6) (with the condition at infinity) was already considered in sec¬ 
tion 


5.1 and can be exactly eliminated using the discrete impedance function (5.1) 
1 f Uj^i — Uj Uj — Uj_i 


h 


-Auj = qj, j = 

- Auq = 0, M-^-i = 0. 


This formally corresponds to a Schur complement. Upon substitution Rn{A) « Fh{A) 
we arrive at the approximate problem 


l/M]Vl-M" <-<-1 


+ l “J 

h \ h 


0.5h 


-i?„(A)-i<- 


J 3- 

h 

Uq — 


- Au^ = q, , 


- AuJ^ = 0 , 


j = -£,...,-l, 
m”, 1 = 0. 


Hence 


|Mf-M,|| =0(||i?„(A)-U,(A)||), 


since all the involved linear systems are well posed uniformly in n. 

Performing similar manipulations with the approximate problem in reverse order, 


we obtain the equivalent system (5.7)-(5.8) 




0.5/1 


-b - 


U^ — 


AUj — qj, j — 1, 


- Au^ = 0, 


(5.7) 


^Problem (5.3i—(5.4i can be viewed as the second-order discretization of ^ 2 '*^ ~ 

= 0, wH! , _. = 0 for some regular enough q supported on [—h{i -|- 1), 0]. As the infi- 


I x= — h(£+l) 


I 2I = -l-00 


nite exterior problem ||5.6|| approximates with second-order ac cura cy the same equation on [0, +co) 
with conditions = —b and = O 5 the relation (5.5i approximates with second order 

the same equation restricted to 1 ), 0 ] with conditions = 0 and = —b. 
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1 


+ l '^0 

hj+i 


b - AuS = 0 , 


S - 


- Au^ = 0 , 


(5.8) 


j = 1, 


by introducing b and fictitious variables m" with positive subindices which, unlike 
their negative counterparts, do not approximate corresponding components of uix). 


Finally, eliminating b we can merge the systems (5.7)-(5.8) into a single recursion 


1 


’^j+i 

hj+i 




Qjj 


j = n-1, 

= 0, rtF = 0, 


^-e 


with the convention that hj := h for j < 0, hj := h for j < 0, 


hj := hj for j > 0, and hg := ho + h/2 (see also Figure 5.1). This finite difference 


hj := 


hj for j > 0, 


scheme is easy to implement by simply modifying the n trailing primal and dual grid 
steps in a given finite difference scheme with step size h. We reiterate that this scheme 


converges exponentially with error 0{\\Rn{A) — Fh{A)\\) to the solution of (5.3|-(5.4| 
in the interior domain, i.e., for the nonpositive subindices. 

The above derivation can easily be extended to variable operators A = Aj in the 
interior domain and tensor-product PML discretizations. This will be illustrated by 
a numerical example in section [6.2| 


Dual steps 


“T 


h/2 Hq 

-1- 




Primal steps h h hi /12 ^3 


Figure 5.1. Schematic view of a finite difference grid appended with an absorbing boundary 
layer generated from quantities in the continued fraction ( |3.3| l . The example shown here is for the 
case n = 3. The gray-shaded region corresponds to the appended absorbing boundary layer, and the 
grid steps ho, hi, ..., hn-i and hi, h 2 ,... ,hn in this layer are generally complex. 

6. Numerical experiments. 

6.1. Waveguide example. To test the accuracy of our absorbing boundary 
layer, we consider the inhomogeneous Helmholtz equation 

Au{x, y) + k'^u{x, y) = f{x, y) 

on a rectangular domain H = [0, L] x [0, H] of length L and height H. We prescribe 
homogeneous Dirichlet conditions at the upper and lower boundaries in y. The source 
term is set to 


/(x, y) = 10 • S{x — 51l7r/512) • S{y — 507r/512), 


with the Dirac delta function 6{-). 
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Our aim is to verify that our absorbing boundary layer models the correct physical 
behavior. To this end we solve the above Helmholtz equation on two rectangular 
domains with fixed height H = tt and different lengths L = n and L = 27r, respectively. 
See also Figure [6T] (left and right, respectively). The wave number is chosen as k = 50. 
The problem is discretized by central hnite differences with step size h = 7r/512 in 
both coordinate directions. The eigenvalues of the resulting tridiagonal matrix A, 
corresponding to the operator —d’^jdy^ — on [0,7r] with homogeneous Dirichlet 
boundary conditions, are explicitly known and eigenvalue inclusion intervals are 

[oi, 6i] U [ 02 , & 2 ] = [—2.50e3, —1.95el] U [7.98el, 1.04e5]. 


We extend the interior finite difference grid by our absorbing boundary layer with 
n = to/ 2 additional grid points to the left of a; = 0 and to the right of a: = L, with 
the near-optimal grid steps computed from a rational interpolant Rn{z) of Fh{z) as 
explained in section |5.1| The physical domain can hence be thought of as an infinite 
strip parallel to the a;-axis. We therefore expect the solutions of both problems (with 
L = TT and L = 27r) to coincide when they are restricted to [0,7r] x [0, tt ]. Visually, 
this is indeed the case, as one can see in Figure 6.1 (where n = 10). Note how the 
amplitude of the solution is damped very quickly inside the absorbing boundary layer. 

To quantify the accuracy of our absorbing boundary layer numerically, we plot in 
Figure [6/^ the relative uniform norm of the difference of the two numerical solutions 
ui(x,y) and U 2 {x,y) restricted to [0,7r] x [0,7r], i.e.. 


err = max \ui{x,y) — U 2 {x,y)\/ max |Mi(a;,y)|. 

0<ai,y<7r / 0<a;,i/<7r 


( 6 . 1 ) 


Indeed, this figure reveals exponential convergence with the rate p given in Theo¬ 
rem [2^ In this example, the expected rate \s p ~ 0.57 and this is indicated by the 
slope of the dashed line in Figure [6(^ 

We would like to mention that absorbing boundary layers usually require some 
physical separation from the support of the right-hand side (the source term) [32| . 
However, thanks to the efficient absorption of evanescent and propagative modes even 
on spectral subintervals with extreme interval ratios, we are able to place our Dirac 
source extremely close to the PML boundary (only one grid point away, see the right 
of Figure 6.11 without deteriorating convergence (see Figure [6(^. 


6.2. PML in multiple coordinate directions. In this experiment we demon¬ 
strate how our perfectly matched layer can be used to mimic domains which are 
unbounded in several coordinate directions, and where there is a nonconstant wave 
speed. To this end consider 


c{x, y)'^Au{x, y) + k'^u{x, y) = f{x, y) 


on a square domain Hi = [0,1]^. The wave speed c(x,y) varies as indicated in 
Figure |l.l[ with c = 1 in the gray region (background material), c = l/-\/2 i n th e light 
gray layer, and c = \/2 in the dark inclusion at the bottom (see also Figure The 
wave number is chosen as k = 120 and the source term is set to 


fix, y) = 5ix - 120/400) • 5iy - 280/400). 

The domain Hi is discretized by central finite differences with step size h = 1/400 
in both coordinate directions. We aim to append absorbing boundary layers with 
n € {7,9,11,13} grid points at each of the four edges of Hi. 
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amplitude amplitude 




Figure 6.1. Amplitude (top) and phase (bottom) of the solution to the waveguide problem in 
section \6. 1\ on two rectangular domains (left/right) which differ in their length. The left domain is 
of length L = 27r in the x-direction, whereas the right domain is of length L = tt. Both domains have 
been appended with absorbing layers at the left and right boundaries. As the absorbing boundary 
layers serve the purpose of extending the physical domain towards infinity, both solutions are expected 
to coincide on the restriction to x £ [0, tt]. In these pictures we have chosen m = 20, so there are 
n = 10 points appended to the left and right boundaries. The step size in the interior domain is 
h = 77)612 in both coordinate directions. 



Figure 6.2. Exponential convergence of the accuracy of the absorbing boundary layers for the 
waveguide problem in section \6.1\ with varying Zolotarev parameter ni £ {8,12,..., 36} (twice the 
number of grid points in each absorbing boundary layer). The expected convergence rate p ~ 0.57 
by Theorem\2. is indicated by the dashed line. 


For constructing the absorbing layers in the j/-direction (below y = 0 and above 
?/ = 1 ) we need inclusion intervals for the negative and positive eigenvalues of Lx — k'^I, 
where is the discretization of —c{x)'^d^jdx'^ on [0,1] with homogeneous Dirichlet 
boundary and c{x) = 1 for y G {0,1}. (Note that c{x, y) varies only tangentially along 
the boundaries of fli, so for y G {0,1} we can indeed write c{x,y) = c{x).) Possible 
inclusion intervals for the eigenvalues are 


[ai,bi] U [ 02 , 62 ] = [-1.44e4,-2.53e2] U [4.94e2,6.26e5]. (6.2) 

For constructing the absorbing layers in the a:-direction (to the left of x = 0 and 
to the right of x = 1 ) we need inclusion intervals ( 6 . 2 ) for the negative and positive 
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eigenvalues of Ly — k^I, where Ly is the finite difference discretization of — c{y)^jdy'^ 
on [0,1] with homogeneous Dirichlet boundary and 

c( ) = / -y- 

^ |l, otherwise. 

Possible intervals are 


[ai,bi] U [ 02 , 62 ] = [-1.44e4,-2.42e2] U [4.82e2,6.26e5]. 


(6.3) 


From the union of intervals in (6.2) and (6.3) we can now calculate the grid steps 


of absorbing boundary layers in the y- and a;-directions, and then modify the finite 
difference matrices to and Ly, respectively. As in the previous example, this is 


done by computing a rational interpolant Rn{z) of Fh{z) defined in section [5d] 

However, there is a small subtlety one has to be aware of with the approach 
just described: effectively, the NtD operators are now given as Fh{Lx — fc^J) and 
Fh{Ly — respectively, and the involved matrices are no longer Hermitian. In 

Figure 


6.4 


(left) we show the eigenvalues of £ 3 , and Ly, respectively, and observe that 
these eigenvalues have “lifted off” the real axis into the upper half of the complex 


plane, in agreement with the analysis of 16 for continuous one-dimensional damped 


operators. From Figure [24] we find at least visually that the Zolotarev approximant 
is of a good quality in this region as well, and the accuracy of the resulting absorbing 
boundary layers should still be satisfactory. 

To quantify the accuracy numerically, we solve the same Helmholtz problem on 
a smaller domain Q .2 = [0.1,0.9]^, again appended with absorbing boundary layers of 
n grid points at each of the four edges of H 2 . As the source term f{x, y) is supported 
inside H 2 , we expect coinciding solutions ui{x,y) and U2{x,y) on their restrictions 
to H 2 . In Figure | 6 ^ (right) we have plotted the relative uniform norm of the difference 
of both solutions, i.e., 


err = max \ui{x, y) -U2(,x,y)\ max |mi {x, y)\. 

0.1<rE,y<0.9 / 0.1<ai,y<0.9 


(6.4) 


Again we observe exponential convergence, and the reduction of the measured error 
is in good agreement with (even slightly better than) the rate p — 0.59 expected from 
Theorem 12.31 


7. Summary, generalizations, and open problems. We have presented a 
new approach for the construction of discrete absorbing boundary layers for indefi¬ 
nite Helmholtz problems via complex coordinate transforms. This approach is based 
on the use of near-optimal relative rational interpolants of the inverse square root 
(or a modification thereof) on a negative and a positive real interval. Bounds for 
the approximation error have been derived, and the exponential convergence of the 
approximants has been established theoretically and demonstrated at numerical ex¬ 
amples. Although our focus in this paper was on absorbing boundary conditions for 
indefinite Helmholtz problems, it was recently understood that these conditions also 
constitute good approximations to Schur complements of certain PDF discretization 
matrices, and they became a crucial component of modern Helmholtz precondition¬ 
ers, such as Schwarz domain decomposition 12 23 and the sweeping preconditioner 
in |22j . Preliminary results have shown successful application to a multilevel do¬ 
main decomposition preconditioner, and a related Schlumberger patent application is 
pending. 
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Fi gur e 6.3. Amplitude (left) and phase (right) of the solution to the Helmholtz problem in 
section \6. 2\ on a square domain Qi = [0,1]^ appended with absorbing boundary layers at all boundary 
edges. In these pictures we have chosen the Zolotarev parameter m = 14, so there are n = 7 
grid points appended to the boundaries. The step size in the interior domain is h = 1/400 in 
both coordinate directions. The dashed square in the interior indicates the smaller domain 0,2 = 
[0.1, 0.9]^, on which we solve the same Helmholtz problem for assessing the numerical accuracy of 
our absorbing boundary layers. 



Figure 6.4. Left: Eigenvalues of the matrices Lx and Ly associated with the Helmholtz problem 
in section \6.2\ appended with n = 7 grid points at the boundaries. Right: Exponential convergence of 
the accuracy of the absorbing boundary layers with varying Zolotarev parameter m E {14,18,..., 26} 
(twice the number of grid points in each absorbing boundary layer). The expected convergence rate 
p ~ 0.59 by Theorem\2.3\ is indicated by the dashed line. 


7.1. Time-domain problem. Classical (explicit) finite-difference time-domain 
formulations lead to PMLs that can be represented via grid steps and 7 ^ which are 

dependent on the wave number k as = a^ -t- %, 77 = a, -t- %, where aj , /3j , Sj , (3^ are 

PPM 


real positive parameters 


Our experiments suggest that the steps of our op¬ 


timal PML grids always have positive real parts, and negative or zero imaginary parts; 
see, e.g., the grid in Figure 3.1 So formally, the steps of the frequency-dependent 
PML can be obtained as 7 , = -|- 7 , = 


-I- i^hj ^, where the steps 


hj and hj are obtained for a fixed wave number fco- If the rational approximant for 
k = ko uses symmetric intervals of approximation, the corresponding grid lies on a 
semiaxis rotated by — j with respect to K+ and the introduction oi k ^ k^ is equiva¬ 
lent to the rotation of the grid and the spectral measure respectively on C \ IR+ and 
C\M_. Therefore such grids retain the exponential convergence for z G K_. However, 
it is not clear if the exponential convergence holds for nonsymmetric intervals, neither 
is known if this convergence holds for the approximation of the discrete impedance 
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function from section Or. 11 

Alternatively, the time domain solution can be represented via stability-correcting 
functions of the discretized operator with the PML obtained for a fixed wave number; 

It can then be efficiently computed in the time domain via Krylov subspace 


see 
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projection. 

7.2. Maxwell and elasticity systems. Important hyperbolic systems, such as 
isotropic Maxwell’s and linear elasticity systems are usually approximated via stag¬ 
gered finite-difference (finite volume) schemes 53 55 . Staggered schemes for mul¬ 


tidimensional problems can be constructed via tensor products of one-dimensional 
staggered schemes. Thus the one-dimensional staggered grids developed in this pa¬ 
per can be automatically implemented in such framework, similar to what was done 


15 39 for PMLs based on the single interval rational approximants. Work 


also provides error estimates for the propagative modes in isotropic elasticity systems, 
showing that for hyperbolic systems one may need rational approximants on slightly 
larger spectral intervals compared to the scalar wave equation. 

7.3. Adaptive grids. The uniform approximation approach requires bounds 
for the smallest/largest negative and positive eigenvalues, which can be rather loose 
due to the weak dependence of the convergence rate of the Zolotarev approximants 
on the interval ratios. The external bounds of the intervals can thus be estimated 
roughly. Still, the numerical estimation of the internal bounds can be rather difficult, 
and accidentally at least one eigenvalue may be very close to the origin, in which case 
even an optimal approximant may require significant order for a satisfactory accuracy. 
To circumvent this problem, it would be interesting to derive a parameter-free near- 
optimal rational approximant of which takes into account the discrete nature 

of the spectrum of A and the spectral weights of the vector b. Promising steps have 


been made by using adaptive rational Krylov algorithms 11 30 31 for this purpose. 


7.4. Variable coefficients. As explained in this paper, variable PDE coeffi¬ 
cients in the tangential direction can be straightforwardly incorporated into the PML 
by modifying A. Moreover, according to preliminary findings, our PML approach 
may be generalized for coefficients varying in the normal direction. Let us replace the 


first equation of (1.2) by 


-I- c{x)u{x) = Au{x), 

ox^ 


with compactly supported coefficient c{x) G Loo- Then grid steps hj,hj, j = 1,..., n 
can be obtained via rational approximation of the “fine grid” finite-difference NtD 
map , involving the finite-difference system 


1 / Uj + i - Uj 

h\ h 



+ c{jh)uj = zuj, 


j = 0,1,..., 00 


with boundary conditions (5.4 ). F or constant c such an approach yields being the 
same as F^ defined in section 5.1 Our experiments with discrete PMLs for variable 


coefficients exhibited exponential convergence albeit at a slower rate than for the 
constant coefficient problem. 


7.5. Connection to inverse problems. Constructing PMLs can be viewed as 
finding equivalent media matching the NtD maps, and this is reminiscent to what 
is done in inverse problems of electrical impedance tomography (FIT). In fact, the 
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conversion of rational approximations to finite-difference schemes (and its planar gen¬ 
eralization) was the basis for the solution of EIT inverse problems via resistor network 
approximations [^. 

Finally, we would like to point out that cloaking problems (which are popular in 
the inverse problems community) are closely related with the construction of PMLs, 
because the latter can be viewed as cloaking of the point at infinity. Cloaking problems 
can also be formulated via complex coordinate transforms 38 and lead to approxi¬ 


mation problems of NtD maps. Although the involved Stieltjes impedance function 
F{z) is typically different in these applications, techniques similar to those presented 
in this paper may still be applicable. 
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Appendix A. Zolotarev approximation and proof of Theorem 

The solution of the Zolotarev problem (2.4) can be computed as 


i=i 


(c.d).. 


(c.d) 


d ■ dn 


(2m - 2j + 1)K(<5') 
- ^ 


(A.l) 


where 


S = c/d, S' = Vl-(52, 


K(<5) = 

Jo 


0 


dt 


is the complete elliptic integral of the first kincQ and where the Jacobian elliptic 
function dn is defined via another such function, sn, by the relations 


dn(w, k) = -y/f — sn(tt, k), ^ = sn(M; k), 


dt 


u = 




In order to prove near-optimality results, we first need to study th e qua ntity 


in (2.4) carefully. Evidently, < 1. Upper and lower bounds for (2.4) were given 

40) as 


exp 


< < 2 


2 — 


< zexp 


7rK(/r') 
-'1 

4K(/r) 


(A.2) 


with 




^T he definition of K{5) is not consistent in the literature. We stick to the definition used 
in |44| Ch. VI]. In Matlab one would type ellipke (delta'‘2) to obtain the value K(5). 
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Hence the Cauchy-Hadamard convergence rate can be computed as 


m=expf'-!^V 


4K(m) J 


(A.3) 


Recalling the equalities ( |2.6[ ) and ( |2.8[ ), let us define the sets 

K = F{K), Ki=F{Ki), k2=F{K2), 

and consider the following auxiliary problem: find a (complex) monic polynomial 
of degree m being the minimizer of 


mm max 
LfePm sGif 


H{s) 


H{-s) 


(A.4) 


We now construct an approximate solution of this problem and show that the approx¬ 
imate solution gives the maximum in (A.4) which yields the best possible functional 
value up to a moderate multiplier. 

Accounting, as it was done in 17 Section 2], that 

= 1 if s G .^2 


,v 

^mi 


^mi 

'^^\^s) 


and 




is) 




i-s) 


= 1 


if s G Ki 


because these polynomials have real coefficients, the polynomial defined in (2.5) 
satisfies 


max 

seKi 


Hrais) _ 

TT ( \ ~ ^^2 

Hm ( S) 


and 


max 

sGK2 


^mjs) _ p(V-bi,v^=ar) 
Hmi-S) 


Lemma A.l. The polynomial Hm defined in (2.5) satisfies the inequality 
Hm{s) 


max 

sGA 


Hmi-s) 


/ „ r -1/2 -l/2'l n 

< 2max|pi ,p2 jp 


(A.5) 


(A.6) 


(A.7) 


with the numbers pi,p 2 cmd p defined in (2.6) and (2.8), provided that mi,m 2 are 
chosen according to l2. 


On the other hand, for any complex polynomial H G Vm we have 

His) 


Hi-s) 


max 

seK 


> P^- 


(A.8) 
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Proof. Let Hm be defined as in (2.51 and conditions (2.71 be satisfied. Accounting 


for (|A.5|), (|A.6|) and (|A.2|), we obtain 

< 2max{p™\p™"} = 2p'"max{pi, 


max 

s£K 


H^{s) 


Hra{-S) 


which gives assertion dAJl) 

To prove assertion ( |A.8[ ) , we consid er t h e th ird Zolota rev problem in the complex 
plane for the condenser (AT,—AT) (see 25 , 54 § 8.7] or 47, § VIII.3]). Due to the 


symmetry of the condenser, the two measures forming the (unique) equilibrium pair 
for (A', — AT) are symmetric to each other in the evident sense. Thus, one can choose an 
(in the Cauchy-Hadamard sense) optimal sequence of type (m, m) rational functions 
of the form H{s)/H{—s), deg{H) = m > 1, such that the roots Sj (1 < j < m) of 
each polynomial H belong to AT. Define 

77(i)(s)= (s-Sj), deg(i7(^^) = mi, 

l<j<m 

Sj 

H^‘^\s)= [s-Sj), deg(i7^^^) = m 2 , 

l<j<m 

Sj GK 2 

mi + 1712 = 'rn. 


By virtue of (A.2) and the location of the roots we have 


and 


whence 


max 

seKi 


max 

SGA 2 


His) 


Hi-s) 


Hjs) 

Hi-s) 


= max 

seKi 


= max 


SGA 2 


i 7 (i)(-s) 


i7(2)(s) 

77(2) (-s) 


> 


o m 
2pi 


> 


1 I 2mi 

i + Pi 




max 

sGK 


His) 


Hi-s) 


max 

s£-K 


Hjs) 

L^(-s)J 


= max 

s£K 


His) 


Hi-s) 


1 I 

1 + P 2 


> max 


2m2 ’ 


2p” 


1 I 2mi ’ -I I 2m2 

i + Pi 1 + P 2 


> max{pi™\p2'""} ■ 


Since the quantity max , P 2 ^^} under the conditions xi >0, X 2 > 0, X 1 +X 2 = m 
is minimal at 


xi = m- 


log P 2 


log Pi + log P2 ’ 


X 2 = m ■ 


log Pi 


log Pi + log P2 


we obtain 


so 


max 

sGK 


Hjs) 

Hi-s) 


max 

sG-K 


Hjs) 

Hi-s)\ 


-1 


>p" 


as m —>■ 00 , 


liminf max 

m^oo I 


His) 


Hi-, 


max 

sG-K 


H(a 

Hf^) 


1 -1 


l/n 


>p'- 
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It follows in view of 25 Theorem 1, Formula (12)] that the logarithmic capacity of 
our condenser satisfies 


exp (—1/cap (iF, —K)^ > P^- 


Moreover, 25, Theorem 1, Formula (11)] yields for all H,G €Vr, 


max 

seK 


H{s) 


G{s) 


max 

se-K 


G{s) 


H{s) 


>P" 


from which (A.8| follows. 


We are now prepared to conclude the proof of Theorem 2.3 
Proof. To establish (2.10), it suffices to note that 


max 

zGK 


Rn{z) 

F{z) 


= max 

sGK 


sPn-i(s^) 

Qnis^) 


□ 


and to apply (A.7) from Lemma A.l condition (2.9), and a consequence of (2.3) for 
finding 


sPn-ijs'^) 

Qn{s'^) 


2 

H^(s) 




1 , 

" ' H^{-s 



To justify (2.11), set z = s'^ and assume that for some pair {P, Q) and R = P/Q 
we have the inequality 


max 

zGK 


F{z) 1 + p'^ 


Define FI by means of (2.1) and rewrite the equality (2.3) in the form 


H{s) 

H{-s) 


RG) _ 1 

F(z) ^ 


RG) 

F(z) 



+ 2 


We readily derive 


max 

sGA 


H{s) 

H{-s) 


< 


1+P" 

_ 2p^ 

l+p" 


= P 


which contradicts (A.8) and thereby proves the assertion (2.11). 


□ 
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